Spotify CS Forecast Exercise

The Analysis

Using the data set provided and outlined below, we’d like you to build a forecast for customer service email volume at the daily level until December 31, 2021. The purpose of the forecast will be used to determine our daily customer service staffing needs based on the email volume of the day.

Load Libraries

library(xgboost)
library(tidymodels)
library(modeltime)
library(tidyverse)
library(lubridate)
library(timetk)
library(skimr) 
library(readxl)
library(scales)

The Data

CS_Forecast_Exercise_Data <- read_excel("data/CS Forecast Exercise Data.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"))

# check date format
# str(CS_Forecast_Exercise_Data)

# change to Date
CS_Forecast_Exercise_Data <- CS_Forecast_Exercise_Data %>% mutate(date = as.Date(date))

skim(CS_Forecast_Exercise_Data)

CS_Forecast_Exercise_Data %>%
  select(date, email, subs, subs_standard, subs_student, subs_family) %>% 
  filter(date <= '2021-05-31') %>% 
  pivot_longer(!date, names_to="type", values_to="count") %>% 
  ggplot(aes(x = date, y = count)) +
  geom_line() + 
  facet_wrap(~ type, scales = "free_y") +
  scale_y_continuous(label=comma) +
  theme_minimal()

CS_Forecast_Exercise_Data %>% 
  select(date, email) %>% 
  filter(date <= '2021-05-31') %>% 
  plot_time_series(date, email)

CS_Forecast_Exercise_Data %>% 
  select(date, email) %>% 
  filter(date <= '2021-05-31') %>% 
  plot_anomaly_diagnostics(date, email)

CS_Forecast_Exercise_Data %>% 
  select(date, email) %>% 
  filter(date <= '2021-05-31') %>% 
  mutate(email = ts_clean_vec(email)) %>% 
  plot_time_series(date, email)

CS_Forecast_Exercise_Data %>% 
  select(date, email) %>% 
  filter(date <= '2021-05-31') %>% 
  mutate(email = ts_clean_vec(email)) %>% 
  plot_seasonal_diagnostics(date, email, .interactive = FALSE) +
  theme_minimal()

Forecast

Next we’ll set up a forecast framework where various models will be trained and evaluated.

First well split our data into training and testing portions. This will allow us to evaluate our models.

# create forecast series
cs_ts <- CS_Forecast_Exercise_Data %>% 
  select(date, email) %>% 
  filter(date <= '2021-05-31') %>% 
  mutate(email = ts_clean_vec(email))

# Split Data 80/20
splits <- initial_time_split(cs_ts, prop = 0.8)

Next we’ll create multiple models to evaluate

# Model 1: auto_arima ----
model_fit_arima_no_boost <- arima_reg() %>%
    set_engine(engine = "auto_arima") %>%
    fit(email ~ date, data = training(splits))

# Model 2: arima_boost ----
model_fit_arima_boosted <- arima_boost(
    min_n = 2,
    learn_rate = 0.015
) %>%
    set_engine(engine = "auto_arima_xgboost") %>%
    fit(email ~ date + as.numeric(date) + day(date),
        data = training(splits))

# Model 3: ets ----
model_fit_ets <- exp_smoothing() %>%
    set_engine(engine = "ets") %>%
    fit(email ~ date, data = training(splits))

# Model 4: prophet ----
model_fit_prophet <- prophet_reg() %>%
    set_engine(engine = "prophet", weekly.seasonality=TRUE) %>%
    fit(email ~ date, data = training(splits))

# Model 5: lm ----
model_fit_lm <- linear_reg() %>%
    set_engine("lm") %>%
    fit(email ~ as.numeric(date) + wday(date),
        data = training(splits))

Now we’ll add these to a model structure that allows us to more easily handle model evaluation

models_tbl <- modeltime_table(
    model_fit_arima_no_boost,
    model_fit_arima_boosted,
    model_fit_ets,
    model_fit_prophet,
    model_fit_lm
)

models_tbl

# Modeltime Table
# A tibble: 5 x 3
  .model_id .model   .model_desc                             
      <int> <list>   <chr>                                   
1         1 <fit[+]> ARIMA(4,0,0)(0,1,2)[7]                  
2         2 <fit[+]> ARIMA(1,0,3)(0,1,1)[7] W/ XGBOOST ERRORS
3         3 <fit[+]> ETS(M,N,M)                              
4         4 <fit[+]> PROPHET                                 
5         5 <fit[+]> LM                                      

Now we’ll evaluate the models on test data and add this data to our model table structure.

calibration_tbl <- models_tbl %>%
    modeltime_calibrate(new_data = testing(splits))

calibration_tbl

# Modeltime Table
# A tibble: 5 x 5
  .model_id .model   .model_desc                         .type .calibration_data
      <int> <list>   <chr>                               <chr> <list>           
1         1 <fit[+]> ARIMA(4,0,0)(0,1,2)[7]              Test  <tibble [177 × 4…
2         2 <fit[+]> ARIMA(1,0,3)(0,1,1)[7] W/ XGBOOST … Test  <tibble [177 × 4…
3         3 <fit[+]> ETS(M,N,M)                          Test  <tibble [177 × 4…
4         4 <fit[+]> PROPHET                             Test  <tibble [177 × 4…
5         5 <fit[+]> LM                                  Test  <tibble [177 × 4…

Now we can use this framework to visually inspect our models.

calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(splits),
        actual_data = cs_ts
    ) %>%
    plot_modeltime_forecast(
      .legend_max_width = 25, # For mobile screens
      .interactive      = TRUE
    )

All models have wide confidence intervals and don’t seem to capture the decreasing trend beginning approximately April 2021. Both ARIMA models capture the weekly patterns well and deviate the least in the period after April 2021.

Now we should also take a look at our error metrics in tabular form. The boosted ARIMA model has the lowest mean absolute error and highest r squared.

calibration_tbl %>%
    modeltime_accuracy() %>%
    table_modeltime_accuracy(
        .interactive = FALSE
    )

Now we’ll refit these models to the entire data set and forecast forward 7 months to the end of 2021. All models except for our linear model provide reasonable step ahead forecasts for daily email volume to the end of 2021.

refit_tbl <- calibration_tbl %>%
    modeltime_refit(data = cs_ts)

refit_tbl %>%
    modeltime_forecast(h = "7 months", actual_data = cs_ts) %>%
    plot_modeltime_forecast(
      .legend_max_width = 25, # For mobile screens
      .interactive      = TRUE
    )

Conclusion / Next Steps

• All models except for our linear models capture the weekly pattern of email volume well and provide reasonable forecasts
• In testing each model the ARIMA & ARIMA Boost models performed the best
• Choose a model considering both accuracy, interpret-ability, and computational cost
• Develop further studies(s) based on model(s) chosen
• Develop computational workflow based on model(s) chosen

Epilogue

As an example, we’ll assume that a consensus has been made to explore the regression model further. To do this the model is isolated and further developed by adding features such as nominal day of week and holiday. From this we can view our model coefficients for a preliminary confirmation of seasonal effects, as well as those of holidays.

# recipe
lm_rec <- 
  recipe(email ~ ., data = training(splits)) %>% 
  step_date(date, features = c("dow", "month")) %>%               
  step_holiday(date, holidays = timeDate::listHolidays("US")) %>% 
  step_rm(date) %>% 
  step_dummy(all_nominal(), -all_outcomes())

# model specs
lm_spec <- linear_reg() %>%
  set_engine(engine = "lm")

# lm wflow
lm_wflow <- 
  workflow() %>% 
  add_model(lm_spec) %>% 
  add_recipe(lm_rec)

# lm fit
lm_fit <- 
  lm_wflow %>% 
  fit(data = training(splits))

# view fit coefs
lm_fit %>% 
  pull_workflow_fit() %>% 
  tidy(conf.int = TRUE) %>% 
  ggplot(aes(term, estimate)) +
  geom_point() +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  coord_flip() +
  labs(title = "Coefficients of a linear regression model")

We can also check to see if the model has improved with the added features. Using our previous model as comparison we can see that our rmse has improved from a previous value of approximately 2130 down to about 1750. This improvement is along the lines of the best time series model in terms of rmse and r2 measures. We now have two viable choices for forecasting, an ARIMA variant as well as a more parsimonious linear model.

# refit on all
last_lm <- 
  lm_wflow %>% 
  last_fit(splits)

# view metrics
last_lm %>% 
  collect_metrics() %>% 
  select(-.config, -.estimator) %>% knitr::kable()